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The ground-state energy of the Hubbard model on a Bethe lattice with infinite connectivity at half filling is 
calculated for the insulating phase. Using Kohn's transformation to derive an effective Hamiltonian for the 
strong-coupling limit, the resulting class of diagrams is determined. We develop an algorithm for an algebraic 
evaluation of the contributions of high-order terms and check it by applying it to the Falicov-Kimball model that 
is exactly solvable. For the Hubbard model, the ground-state energy is exactly calculated up to order t l2 /U 11 . 
The results of the strong-coupling expansion deviate from numerical calculations as quantum Monte Carlo (or 
density-matrix renormalization-group) by less than 0.13% (0.32% respectively) for U > 4.76. 

PACS numbers: 71.30.+h, 71.10.Fd, 71.27.-Ha, 02.10.Ox 



INTRODUCTION 

THE Hubbard model 1 captures the essential elements of the 
complex behavior of strongly-correlated fermionic sys- 
tems with short-range repulsive interaction. Particularly in- 
teresting is the exploration of the transition from the param- 
agnetic metallic phase to a paramagnetic Mott-Hubbard in- 
sulator in the limit of infinite dimensions 2 " 6 ; there, the inter- 
acting lattice problem can be mapped onto effective single- 
impurity models and solved within the framework of the dy- 
namical mean-field theory (DMFT). Since this phenomenon 
was discussed controversly 7-9 , high accuracy in determining 
the ground-state energy and double occupancy per lattice site 
near the transition region is necessary for resolving doubts as 
to the nature of the transition, for minimizing quantitative un- 
certainties in the phase diagram and establishing an essential 
benchmark for other, in particular numerical methods. There 
were very many attempts to study the model in the strong- 
coupling limit (cf., e. g., refs. 10-12). However, it appears to 
be rather difficult to go beyond the lowest orders. Therefore, 
we developed a computer-algebraic approach. 

In this work we present a detailed description of a "divide- 
and-conquer" algorithm used for an exact calculation of all co- 
efficients in the asymptotic expansion of the ground-state en- 
ergy of the Mott insulator including order t l2 /U n . Results of 
such an algorithm up to f 10 /t/ 9 were already quoted in ref. 13, 
where an extrapolation scheme to infinite order (extrapolated 
Perturbation Theory, "ePT") was introduced. This showed ex- 
cellent agreement with a quantum Monte Carlo (QMC) tech- 
nique, improved the state of the art by 1-2 orders of mag- 
nitude and lead to a well controlled evaluation of the criti- 
cal exponent. Quite recently, our method was applied to the 
Bose-Hubbard model 14,15 . 

The outline of this paper is as follows: In sec. I, we 



show how the effective Hamiltonian is derived following 
Kohn's 16 and Kato's 16 and Takahashi's 17 treatment of the 
strong-coupling limit. Then, in sec. II, the class of diagrams 
defined by the resulting effective Hamiltonian is discussed for 
the Bethe lattice with infinite connectivity, and the algorithm 
for the evaluation of electronic transfer processes on it is de- 
scribed. The concept of this algorithm is ideally suited for 
parallelization that will be done in further work. The results 
are first given for the Falicov-Kimball model that is exactly 
solvable and serves as a test of our treatment (sec. III. A). 
Our main result is given in eqs. (19) and (20) in sec. III.B. 
We apply our method "ePT" (see ref. 13) and compare our 
results with results from DMFT-QMC and DMFT-DDMRG 
(Dynamical Density-Matrix Renormalization Group) 9 for the 
ground-state insulating phase of the Hubbard model. Finally, 
flowcharts that present the essential parts of the algorithm are 
given in the Appendix. 



I. PERTURBATION EXPANSION FOR THE 
STRONG-COUPLING LIMIT 

We investigate spin-1 /2 electrons on a lattice represented by 
the Hubbard model 



H = T + UD, 



(1) 
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where T = — tY,(iJ),a c la c ia * s tne kinetic energy operator de- 
scribing electron hops between nearest neighbour sites i and j 
with the transfer amplitude t, UD = U £j ntt-ny. is the interac- 
tion part including only local contributions n- ia = c' ia c- la . cL 
and c ia are the creation and annihilation operators for elec- 
trons with spin a =|,t on site i- 

In the following, we sketch the calculation of an effective 
Hamiltonian in a strong coupling expansion, as was developed 
in ref. 18 and ref. 17. There it is shown how this expansion 
in 1/(7 is done systematically. The aim is the transformation 
to new particles with an effective Hamiltonian that does not 
change the number of doubly occupied sites. The considera- 
tions are valid for any lattice in any dimension. The operator 



2 



for the kinetic energy T in eq. (1) couples states with a differ- 
ent number of doubly occupied sites. In deriving this effec- 
tive Hamiltonian, a decoupling can be achieved by introduc- 
ing suitable linear combinations. Rotating to such a new basis 
is performed by a unitary transformation = e s developed 
by Kohn 16 for the strong-coupling limit. This transformation 
introduces new particles created by c T so that 



Because 



and therewith 



H(c) = e s ^H(c)e- s ^ = H(c). 



(2) 



(3) 



The generator S is constructed in such a way that the hopping 
of the new particles does not change the effective number of 
doubly occupied sites for the new particles (c), 



[H(c),D(c)}=0. 



(4) 



This requires S (and therefore H(c)) to be an operator series 
in l/U 



!=1 



(5) 



and the unitarity of the transformation implies S 1 " = —S. Ob- 
viously, the wavefunctions can be expressed in terms of the 
new particles, and it follows for the eigenenergies 

(y m {c)\H{c)\y m {c)) = (y m (c)\H(c)\y m (c)) =E m . (6) 

The ground state y/b(c) °f H(c) at half band filling will be 
determined at the end of this section. 

The low orders in l/U are conveniantly found by substitut- 
ing the expansion (5) in (3); to second order in 1 /U: 

H(c)=T (c) + UD(c) + 1 [Si (c) , H(c)} + 

+ ^ [S 2 (c),H(c)} + [Si (c), ^ (c),H(S)}} + ■■■■ (7) 

From the condition (4), the coefficients S„ (c) are determined 
order by order as shown now. The kinetic energy operator 
can be separated in three parts, each of which increases or 
decreases the number of double occupancies by one, or leaves 
it unchanged, 



T(c)=T u + T_ u + T , 



(8) 



where 



Tu = ~t £ "i-a( l -«j,-<r)4 



ict\JCT 

(ij),<T 
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To = ~t 0-- f k,-a-ni- a + 2ni-ani-a)cla d la- 



[TuMc) 



-T v and [T-u,D(c)]=T-u, 



(9) 



(4) is fulfilled when the cross terms Tu and T-u are cancelled 
by [S(c),UD(c)] in the lowest order in U~ l . This is achieved 
by choosing 



S 1 {c) = T u -T- U . 



(10) 



Inserting in (7) and demanding (4), one obtains the condition 
for the next order, i. e., S2 that leads to 



S 2 {c) = [Tu + T-u,To\. 



(11) 



Following this procedure, one determines S(c) order by or- 
der. Since S does not create or annihilate bare particles, the 
vacuum state is equal for both old and new particles. 

The lowest order terms of the resulting 1 /{/-expansion of 
the Hamiltonian H (c) = UD(c) + YZ=o U~% are (here and in 
the following, T = T(c)) 



j 

h=Z p j TS j Tp j> 

j 

~h 2 = '£PjTS)TS)TP j - 



(12a) 
(12b) 



' 1 (12c) 
- ^L( P J Ts2 j TP j TP i+ P i TP J TS ^ TP j) - 



where Pj projects onto the subspace with j > double occu 
pancies and S k j is defined as (Dj = j) 

Pi 



¥J [Uj 



k>0. 



(13) 



Calculation of the next orders in this way needs an increas- 
ing computational effort. Therefore, a computer program has 
been developed that evaluates the general formula of Kato, 
ref. 18 (cf. eq. 22), up to a given order. 

H(c), the first terms of which are given by (12) is the de- 
sired effective Hamiltonian. It is valid in any subspace with 
fixed number of j double occupancies (of particles corre- 
sponding to c^). For large U, the ground state \jfo(c) of H (c) 
must have the lowest value of UD(c), and because H(c) leaves 
the number of doubly occupied sites unchanged, see (4), this 
state does not contain any double occupancies at half band 
filling, so we put j = in the following. So far, the consid- 
erations are valid for arbitrary dimension. Now, we focus on 
the case of half-filling in infinite dimensions. Then, all global 
singlets are ground states of H (c), cf. ref. 19. Therefore, each 
lattice site is equally likely to be occupied by an electron with 
spin f or I, irrespective of the spin on any other lattice site. 
This enables us to perform the ground-state expectation val- 
ues in (\ffo(c)\h n \\jfo(c)) in the case of a half-filled band in a 
second computer program. 

Both these computer programs are of algebraic nature 
(work with integers only) and thus give exact results for any 
given order in l/U. 



II. GRAPHS AND ALGORITHM 

In this section, we give more details regarding the evaluation 
of H{c). In fixed order i, the operator hj containes all possible 
electron hops resulting from i + 1 applications of T; it con- 
tributes to the energy in the order l/U l . In the following, we 
will deal only with states with zero doubly occupied sites (in 
the new particles), so we drop the index and denote Sq = S k 
and Pq = P. A sequence of electron hops described through 
an operator chain PTS k ' ■ ■ S k 'TP is called process 19 . Only 
processes that restore the initial state contribute to the ground- 
state energy. Consider now a given process and perform the 
sum on lattice sites in the operators T. The individual terms 
in this sum are called "sequences" or "diagrams". Because of 
the linked-cluster theorem, we need to keep only connected 
diagrams. Each of these contains n = 2 + ^ sites, connected 
through i + 1 jumps. Diagrams of odd order in t do not con- 
tribute at half filling for any lattice type. Now, we specialize 
our considerations to the Bethe lattice. There, all closed paths 
are self-re tracing. This can be seen in fig. 1 where the possible 
sequences (of hops) for the low orders i = 1 (n = 2) and i = 3 
(n = 3) are displayed on a Bethe lattice of connectivity 3. In 
the following, we put t = t * / \/z (Z is the number of nearest 
neighbours) and consider the limit Z — > °° for fixed t* = 1. 
This limit implies that in the energy, in each order in U~\ 
the leading order t' +l is taken into account. Thus, diagrams 
with more than two transitions between any two given sites 
are suppressed at least as 1 /Z: every additional connection of 
already doubly connected sites is smaller by 1 /Z compared to 
those with only two jumps between any two sites. Since the 
paths are self-retracing, they can be collapsed into 'Butcher 
trees' 20 as also indicated in the right of fig. 1. The position of 
the first hop in the sequence (the index j of ci a at the rightmost 
T) defines the "root" of the tree, cf. fig. 1 . Because there are 
exactly two hops between neighboring sites (Z — > °°), there is 
a one-to-one correspondence between sequences and Butcher 
trees. The number of Butcher trees built with n vertices, A [n], 
is still moderate for moderate n; it is given, see ref. 21, by the 
following recursive definition 

A M : = —jLMn-j] £ (d k (j)A[d k (j)}) , (14) 

with A[l] := 1 and A [2] := 1, and d k (j) is the k th element of 
the set of divisors of j. In table I, we give the number of trees, 
the number of initial spin configurations, and the number of 
sequences for given order of perturbation theory. 

Illustrating the complexity of the problem, figure 2 shows 
the Butcher trees up to seven vertices (representing the con- 
nected sites) and thus all graphs contributing up to eleventh 
order in the perturbation series. 

In order to illustrate our procedure, we show all intermedi- 
ate states for all sequences for all processes in third order in 
fig. 3. In the main part of the figure, these states are displayed 
on a Bethe lattice with connectivity 3. Note that only three 
sites are affected by the hops, the two in the center and the up- 
per right one on this segment of the lattice. All other sites of 
the lattice are unaffected by the hops and we can restrict our 
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: ®— . (n=2) 




FIG. 1 : Correspondence between sequences of hops on the lattice 
(left) and Butcher trees (right): Shown are the sequences in first (n = 
2) and third order (n = 3) and the related Butcher trees. The first hop 
defines the root of the tree (encirceled). 
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20 


596 


5644880 


11 


7 


48 


2712 


1099056000 



TABLE I: Number of trees with n = (i + 3)/2 sites or order i in per- 
turbation theory, A[n], see text. B[n] is the number of the initial spin 
configurations; C[n] is the number of all hopping sequences, i.e. the 
number of different ways to realize all processes of the given order 
with initial spin configurations. 

considerations to the "Butcher trees" shown in the left part of 
the figure. There, the sites of the first hop, the roots of the 
trees, are encircled. The processes are generated by follow- 
ing parts of the hamiltonian (terms containing the sequence 
. . . PTP . . . vanish at half band filling (hf)): 

hf = PTS 1 TS^TS^TP- PTS 2 TPTS l TP. (15) 

This expression coincides with the fourth order (highest avail- 
able) in ref. 10, eq. (6). The probabilities for the occurrence 
of the processes in the paramagnetic phase, multiplied by the 
prefactors of the processes in h}f (1 and —1 here), are given in 
the right column of figure 3. Their sum yields the contribution 
to the ground-state energy. 

The numerical algorithm to calculate the expectation value 
of the operators hj is based on this diagrammatic approach: 
After constructing all / th order Butcher trees for the lattice, 
all possible sequences on them resulting from different terms 
of the hj are generated through a recursive procedure, within 
which the conditions for the realisation of the electron trans- 
fers, as the fulfillment of the Pauli principle and the consider- 
ation of preceding hopping steps on the branches, are defined. 
The first electron hop starts from the root of the graph (encir- 
cled sites in figures 2 and 1) to a neighbour site. In a single 
loop of the program, possible following jumps are tested by a 
subroutine and executed where applicable. Thereby, the sec- 
ond and last electron hop on a branch has to be performed by 
the same spin species as the first one. This requirement guar- 
antees the restoration of the initial spin configuration. The 
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FIG. 2: Butcher trees up to the order i = 11. The order expressed by the number of sites is i = 2n — 3. The encirceled sites denote the root of 
the trees. 



actual number of double occupancies that enters the operators 
S k , (13), is stored and used for the computation of the factor 
for the given process. The final spin configuration determines 
the factor's sign, (— where P is the number of permuted 
spin pairs. Summation yields the contribution of given order 
to the ground-state energy of the Mott insulator. 

As shown now, this algorithm has been successfully tested 
by computing the ground-state energy of the exactly solvable 
Falicov-Kimball model, a simplified Hubbard model with one 
immobile spin species. 



HI. RESULTS 



A. Falicov-Kimball Model 

For the Hamiltonian of the half-filled Falicov-Kimball model 
we refer to van Dongen's fundamental work 22 . We calculated 
with our procedure the ground-state energy on the Bethe lat- 
tice with infinite connectivity (bandwidth W = 2\/2t*) up to 
G(t xl /U~ n ). Taking t* = 1 as our energy unit in the follow- 
ing, we find that all contributions in the series in 1 /U vanish, 



except the first: 



E l K (u) = -— + 0\- T , 



1 

4U 



1 

C/" 13 



(16) 



Next, we verify 23 our result (16) using the exact solution in 
ref. 22. We start from the expression of the kinetic energy 
eq. 4.7; we denote that by E ' T (U). All ground-state en- 
ergies are given as densities (intensive thermodynamic vari- 
able). We use their spectral representation in order to express 
the Green function in eq. 4.7 in terms of the density of states 
z(e,U), eq. 7.5 in ref. 22. Thus 

E* K - T (U) = -2 JdeJ de'z(e,U)z(e',U)j^. (17) 
o o 

Finally, we show by numerical integration for different 
choices of U between 2 and 10 that 



E FK . T ([/) = --L(] ,,(io 



-16' 



(18) 



and that confirms our result, eq. 16. (Here, 10 16 is the nu- 
merical accuracy.) 
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1 

FIG. 3: The intermediate states for the two processes and all sequences contributing to the ground-state energy in third order in jj. The arrows 
correspond to the application of T . The symbol @ denotes a doubly occupied lattice site; the symbols O an d O denote a hole and a singly 
occupied lattice site, respectively. In the right column, the contribution of each sequence to the ground-state energy is indicated. They sum up 



We have to conclude that all higher order hopping contri- 
butions to the ground state energy cancel. The reason may 
be that only one spin species can hop in the Falicov-Kimball 
model. 



B. Hubbard Model 

The calculation of the ground-state energy of the Hubbard 
model to the 1 1 th order yields 



1 1 1 1 19 1 593 1 

23877 1 4496245 1 ( \ 

128 TP 2048 U 11 + [u^ 



(19) 



Consequently, the double occupancy of the original particles 



j-D(U) =dE(U)/dU is given by 



-D H (U) 



11 3 1 95 1 4151 1 
2U 2 ~ + 2U* + YU^ + 
214893 1 49458695 1 



128 U w 2048 U 12 



32 t/ 8 

1 

7JT4 



(20) 



These perturbation-theoretical (PT) results are shown as solid 
lines in figure 4. The comparisons of the first (second) and 
third (fourth) order PT (dotted/dashed lines) demonstrate a 
fast convergence for the values of U shown. The agreement 
with QMC (circles) and DDMRG (crosses) results extrapo- 
lated to zero temperature is excellent for U > 5 (smaller than 
the line width in fig. 4). As U decreases, devations from 
these (numerical) DMFT data increase noticeably, since re- 
sults of finite order PT rapidly become inaccurate as U ap- 
proaches U C[ , the critical interaction. In the following, we de- 
scribe a method of how to estimate the critical coupling U Cl . 
We assume that the radius of convergence of the l/U ex- 
pansion of the energy coincides with the critical coupling 
U Cl , beyond which the insulating phase becomes stable. We 
perform an extrapolation of the computer-aided high-order 
evaluation to infinite order (ePT 13 ) that exceeds former ac- 
curacy in ref. 13. With E H (U) = Y*7=\ a 2sU l ~ 2s , we have 

U Cl = lim. v ^oo V fl 2i+2/fl2.s- In figure 5, we plot v / ozhV«2s 
against 1 /s. As seen in the figure, the data points are fitted 
by a nearly straight line as a function of 1 Taking a slight 
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FIG. 4: QMC and DDMRG 9 results for the ground-state energy E 
(top) and the double occupancy D (bottom) in comparison with PT, 
see eq. (19) and (20). 



curvature into account in a least-squares fit, 



&2s+2 



2s 



u 2 



(2s) 



2 ' 



(21) 



one finds U Ci =4.76, u t = -16.471257, and u 2 = 5.7147072. 
The critical exponent (for details see ref. 13) defined by 
£ C rit(£/) « (U - U Cl ) T ~ l is obtained with T = 3.46, that gives 
support to our assumption 13 for T = 7/2. 

The ePT estimates for the energy are strongly supported by 
QMC results 13 : E PT is converged within 0(10 5 ) for U > 
6, while the ePT provides an estimate for E with a precision 
of the same order above the stability edge of the insulator. 
These ePT results for E have been reproduced at U = 4.8, 
5.5, and 6 within <^(10~ 6 ) using the self-energy functional 
approach/dynamical impurity approach (SFT/DIA) 24 . 

Other methods based on DDMRG give U Cl ~ 4.45 and T = 
2.5 see ref. 7,9. As seen from figure 4, a high accuracy of data 
is indispensable for a correct analysis of the transition. 



CONCLUSIONS 



Using Kohn's unitary transformation, an l/t/-expansion for 
the Hubbard model was derived up to order (1/t/) 11 at zero 
temperature. The expansion has been formulated in terms 
of diagrammatic rules for the calculation of the ground-state 
energy and the resulting double occupancy. These rules re- 
duce the calculation of finite-order contributions to an alge- 
bra which becomes increasingly complex for higher orders. 




0.30 



FIG. 5: Construction of ePT: PT values for U c (s) = \J ai s +ll a ls 
(squares) are extrapolated to l/s — > using a quadratic least squares 
fit (solid line). Evaluations at smaller 1 /s (circles) define ePT coeffi- 
cients to all orders, cf. ref. 13. 



Any step of the rules is carried out exactly by our computer 
program. Explicit results were obtained for the ground-state 



6 energy up to 11 order in l/U. 



An inspection of the contributions of the diagrams shows 
that there are groups of dominant ones, namely the widespread 
diagrams (first ones in each order in fig. 2), and they are sig- 
nificant in view of the metal-insulator transition. This should 
be analysed quantitatively in the large order limit. Then, even 
an exact determination of, e. g., U C{ might be possible. 
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Appendix: FLOWCHART OF THE ALGORITHM 

The kernel of the program is the recursive procedure hopping 
which is calling the procedure hopping_f wd and vice versa; 
their flowcharts are shown in figure 6. 



^begin) — ► 



hopping 

HI 



— END J 



hopping_f wd 



The main task of these "(non-linear) indirect-recursive" 
procedures is to find all electron hopping processes generated 
by the Hamiltonian, which occur on a set of graphs, as shown 
in the preceding section. The graphs are defined in such a 
way that the set of all possible processes of the given order is 
identical to the set of processes starting from the root of the 
graphs. All possible spin configurations on the graphs are gen- 
erated and stored in arrays. Due to symmetry, it is sufficient 





hopping (hop , branch , doubleocc , hole , down , down , 
spinconf ig , jumphist .true , doublehist) 



hopping Chop .branch, doubleocc+1 , hole .pair , down , 
spinconf ig , jumphist , true , doublehist) 



hopping Chop .branch , doubleocc , up .pair , down, 
spinconf ig , jumphist .true , doublehist) 



hopping Chop , branch , doubleocc , down , pair , up , 
spinconf ig , jumphist , true , doublehist) 



hopping (hop, branch, doubleocc- 1 .down , up , up , 
spinconf ig, jumphist , true , doublehist) 



hopping (hop, branch, doubleocc- 1 , up .down, down, 
spinconf ig, jumphist , true , doublehist) 



^ EHD ^ 
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to occupy the root always by an up-spin leading to 2"~ con- 
figurations, where n is the number of vertices of the graphs. 
For the description of the procedures the following variables 
are used, cf. fig. 6: 

hop: number of the current hop 

branch: number of the branch, on which the hop occurs 

doubleocc: number of double occupancies 

pf , pt, ps: spin states (hole, up, down, pair): pf spin state 
of the site out of which the jump occurs, pt spin state 
of the site to which the jump occurs, ps jumping spin 

spinconf ig: structure describing the current spin configura- 
tion on the graph 

jumphist: table containing the jump history on the branches 

direction: direction of the jump, defined by description of 
the graph: true jump forward, false jump backward 

doublehist: sequence of digits representing the history of 
the number of double occupancies 

branch_count: total number of branches in the graph 

The initial call of the procedure hopping is done with 
the following parameters: hop=0, branch=0, doubleocc=0, 
pf=up, pt=down, ps=up, spinconfig, jumphist, true, 
'0'. The procedure executes jumps on all branches 
in both directions; the possibility of a jump is tested 
through the procedure hopping_fwd. Its call is done 



with hopping_fwd(branch_nr ,hop+l ,true) (forwards) 
or hopping_fwd(branch_nr,hop+l .false) (backwards) 
and it checks if on the current branch a hop has already oc- 
cured in the given direction. If not, depending on the spin state 
of the involved sites, the procedure hopping is called, the ar- 
ray of spin states is refreshed, and a next branch is tested. The 
second jump on a branch has always to be performed by the 
same spin species as the first one. 

The recursion has the property that, in case of exiting the 
procedure when the jump was not possible, the values of 
variables resulting from preceding steps are automatically re- 
stored. 

This construction of the algorithm guarantees that all pos- 
sible variants of electron jumps are tested in accordance with 
the assumptions, and that the final spin configuration equals 
the initial spin configuration; therefore this condition does not 
need not to be tested. After the last step (carrying the number 
2*branch_count), the characteristic factor 

for a process appearing in order l/U l is computed. In (A.l) P 
denotes the number of permuted spin pairs, /V sp i n is the num- 
ber of spin configurations not changed by the process. The 
sum runs over all terms of the hamiltonian which generate 
the process and f m is the related factor obtained from equa- 
tion (7) and calculated by a separate algorithm, dj are the 
numbers of double occupancies and kj are also obtained from 
equation (7). The sum of the factors yields the final result. 
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